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Directional dark matter detectors will be able to record the recoil momentum spectrum of nuclei 
hit by dark matter WIMPs. We show that the recoil momentum spectrum is the Radon transform 
of the WIMP velocity distribution. This allows us to obtain analytic expressions for the recoil 
spectra of a variety of velocity distributions. We comment on the possibility of inverting the recoil 
momentum spectrum and obtaining the three-dimensional WIMP velocity distribution from data. 



. I. INTRODUCTION 

o: 

' The identification of dark matter is one of the major open questions in physics, astrophysics, and cosmology. 
Recent cosmological observations together with constraints from primordial nucleosynthesis point to the presence of 
O I non-baryonic dark matter in the universe. The nature of this non-baryonic dark matter is still unknown. 
D One of the preferred candidates for non-baryonic dark matter is a weakly interacting massive particle (WIMP). 

Substantial efforts have been dedicated to WIMP searches in the last decades 0|. A particularly active area |^] are 
' WIMP direct searches, in which low-background devices are used to search for the nuclear recoil caused by the elastic 
CnJ : scattering of galactic WIMPs with nuclei in the detector In these searches, characteristic signatures of a WIMP 
signal are useful in discriminating a WIMP signal against background. 
' A WIMP signature which was pointed out very early |4| is an annual modulation of the direct detection rate caused 
by the periodic variation of the Earth velocity with respect to the WIMP "sea" while the Earth goes around the Sun. 
1*"^ The typical amplitude of this modulation is 5%. A modulation with these characteristics was observed by the DAM A 
collaboration [|| , but in light of recent results ||, Q , its interpretation as a WIMP signal is currently in question. 
^\ Different, and possibly clearer, WIMP signatures would be beneficial. 

A stronger modulation, with an amplitude that may reach 100%, was pointed out by Spergel in 1988 [^|. Spergel 
noticed that because of the Earth motion around the Sun, the most probable direction of the nuclear recoils changes 
, with time, describing a full circle in a year. In particular this produces a strong forward-backward asymmetry in the 
angular distribution of nuclear recoils. 
CIh' Unfortunately it has been very hard to build WIMP detectors sensitive to the direction of the nuclear recoils. A 
promising development is the DRIFT detector The DRIFT detector consists of a negative ion time projection 
^ chamber, the gas in the chamber serving both as WIMP target and as ionization medium for observing the nuclear 
f~| : recoil tracks. The direction of the nuclear recoil is obtained from the geometry and timing of the image of the recoil 
^ ' track on the chamber end-plates. A 1 m^ prototype has been successfully tested, and a 10 m'^ detector is under 
. ^ , consideration. 

In addition to merely using directionality for background discrimination, what can be learned about WIMP prop- 
H \ erties from the directionality of WIMP detectors? It is obvious that different WIMP velocity distributions give rise 
. - - 1 to different recoil distributions in both energy and recoil direction. Copi, Heo, and Krauss [Q, and then Copi and 
Krauss have examined the possibility of distinguishing various WIMP velocity distributions using a likelihood 
analysis of the resulting recoil spectra, which they generated through a Monte Carlo program. They have concluded 
that a discrimination among common velocity distributions is possible with a reasonable number of detected events. 

Here we want to gain insight into the properties of the nuclear recoil spectra in energy and direction. For this 
purpose, we develop a simple formalism that relates the WIMP velocity distribution to the distribution of recoil 
momenta. We find that the recoil momentum spectrum is the Radon transform of the velocity distribution (see 
Eq. (1^) below). We apply this analytical tool to a series of velocity distributions, and discover for example how the 
recoil momentum spectrum of a stream of WIMPs differs from that of a Maxwellian velocity distribution. With our 
gained insight, we suggest that if a WIMP signal is observed in directional detectors in the future, it may be possible 
to invert the measured recoil momentum spectrum and reconstruct the WIMP velocity distribution from data. 

In Section || we describe the general kinematics of elastic WIMP-nucleus scattering, and in Section [II we obtain 
our main formula for the nuclear recoil momentum spectrum. Sections IV and ^ contain general considerations and 



examples of Radon transforms of velocity distributions. Finally, Section VI discusses the possibility of inverting the 
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recoil momentum spectrum to recover the WIMP velocity distribution. The Appendices contain useful mathematical 
formulas for the computation and inversion of 3-dimcnsional Radon transforms. 



II. WIMP-NUCLEUS ELASTIC SCATTERING 

Consider the elastic collision of a WIMP of mass m with a nucleus of mass M in the detector (see Fig. Let 
the arrival velocity of the WIMP at the detector be v, and neglect the initial velocity of the nucleus. After the 
collision, the WIMP is deflected by an angle 9' to a velocity v', and the nucleus recoils with momentum q and energy 
E = /2M. Let 9 denote the angle between the initial WIMP velocity v and the direction of the nuclear recoil q. 
Energy and momentum conservation impose the following relations: 

1 2 1 /2 , 9^ n\ 
— mv — —mv H , 1 

2 2 2M' ^ ' 
mv' cos 9' — mv — q cos 9, (2) 

mv' siTi9' — qsm9. (3) 

Eliminating 9' by summing the squares of Eqs. and (^), 

,2„./2 _ j-^j^ _ qcos9)'^ + (qsin^)^ = m^w^ - 2mvqcos9 + q^ , (4) 



m V ' 



and using this expression to eliminate v' from Eq. (|l|), gives 

q = 2fivcos9, (5) 

where 

mM 

is the reduced WIMP-nucleus mass. We deduce that the magnitude q of the recoil momentum, and the recoil energy 
E, vary in the range 

< g < w = < ^ < S,„ax = (7) 

iVl 

Eq. (|^) will be exploited in the following Section to express the recoil momentum distribution in a simple mathe- 
matical form. For this purpose, we also need the expression for the WIMP-nuclcus scattering cross section. We write 
the differential WIMP-nucleus scattering cross section as 



dfj <tq 
dq^ q 



— = ^S{q), (8) 



where ctq is the total scattering cross section of the WIMP with a (fictitious) point-like nucleus, and S{q) — |F(g)p 
is a nuclear form factor normalized so that 5'(0) = 1. (Both S{q) and F{q) are confusingly called form factors.) 
Eq. ^ is valid for both spin-dependent and spin-independent WIMP-nucleus interactions, although cto and F{q) have 
different expressions in the two cases. For example, for spin-independent interactions with a nucleus with Z protons 
and A~ Z neutrons, 

a^^^[ZGl + {A-Z)G:]\ (9) 

TT 

where Gf and G" are the scalar four-fermion couplings of the WIMP with pointlike protons and neutrons, respectively 
(see Ref. [Q). If the nucleus can be approximated by a sphere of uniform density, its form factor is 

<d[B\n{qR) - qRcos{qR)f 
^^''> - ' ^ ^ 

where 

R ~ [0.91^1/3 -I- 0.3] X lO^^^cm (11) 

is (an approximation to) the nuclear radius. More realistic expressions for spin-independent form factors, and formulas 
for spin-dependent cross sections, can be found, e.g., in Refs. [|l2l O, O, O. 
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III. RECOIL MOMENTUM SPECTRUM 



Eqs. (|) and (|) can be combined to give the differential recoil spectrum in both energy and direction, i.e. the recoil 
momentum spectrum. We define it as the double differential event rate, in events per unit time per unit detector 
mass, differentiated with respect to the nuclear recoil energy E and the nuclear recoil direction q, 

(12) 



dEdQq 



where dQq denotes an infinitesimal solid angle around the direction q. 

The double differential rate follows from the double differential cross section 



dq^dilg 

first through the change of differentials dq^ — 2MdE, and then through multiplication by the number N of nuclei in 
the detector, division by the detector mass MN, and multiplication by the flux of WlMPs with velocities v in the 
velocity space element d'^v, 

nvf{v)d^v. (14) 

Here n = p/m is the WIMP number density, p is the WIMP mass density, and /(v) is the WIMP velocity distribution 
in the frame of the detector, normalized to unit integral. 

The double differential cross section is obtained as follows. Azimuthal symmetry of the scattering around the WIMP 
arrival direction gives dflq = 2t: dcos 9. The relation between cos 9 and q in Eq. (|), cos 6* = q/2pv, can be imposed 
through a Dirac S function, S{cos9 — q/2pv). Thus 



da 



^ xf a ^ \ craSiq) ^( „ q\ 

~ —7 — COS 9 — TT- V COS 9 

dq^ 27r V / 8np,^v \ 2p J 



dq^dflq dq^ 2n 

This is correctly normalized as can be seen by integration of the expression in the middle over dflq . 
Summarizing, the double differential event rate per unit time per unit detector mass is 

^ f „ j3„ naoSiq) f J „ q 



(15) 



dEd^q 



We write it as 

dR naoS{q) 



dEdilq Anp? ' 

Here 



/(«g,q)- (17) 




(18) 

is the minimum velocity a WIMP must have to impart a recoil momentum q to the nucleus, or equivalently to deposit 
an energy E = g^/2M, as can be seen from Eq. (|^). Moreover, 

/(w, w) = y" 5 (v w - u;) /(v) d^v, (19) 

is the 3-dimensional Radon transform of the velocity distribution function /(v). We note in passing that / has units 
of inverse speed. 

Eq. ( ^t|) is the main result of this paper. It states that, apart from a normalizing factor, the recoil momentum 
spectrum is the Radon transform of the WIMP velocity distribution. The Radon transform is a linear integral 
transform (see Refs. |l^, |l^), which was introduced in two dimensions by Radon in 1917 [|8). The Radon transform 
has been widely studied for its use in solving differential equations, and especially in two-dimensions, for its medical 
applications in computer tomography. Geometrically, /(w, w) is the integral of the function /(v) on a plane orthogonal 
to the direction w at a distance w from the origin. For reference, some mathematical properties of the Radon transform 
are given in the Appendices. 
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As a check of our formalism, we show that integrating our basic equ ation ( |17| ) over recoil directions reproduces the 
usual expression for the recoil energy spectrum dR/dE. Applying Eq. ( [A12| ) in Appendix A to our expression for the 
differential rate, we find 

dR _ naoS{q) f /(v) 3 



- 9 2 / (20) 
This is the usual expression of the recoil energy spectrum (cfr. Eq. (8.3) in Ref. 0). 

IV. COMPUTING THE RECOIL MOMENTUM SPECTRUM 

We have cast the nuclear recoil momentum spectrum in terms of a Radon transform. Now we can take advantage of 
the properties of Radon transforms, some of which are listed in the Appendices, to compute recoil momentum spectra 
analytically. In this Section we give some general considerations, and in the next Section we give explicit examples of 
analytic recoil momentum spectra. 

A. Isotropic distributions 

When the WIMP velocity distribution is isotropic, /(v) — f{v), the recoil spectrum is also isotropic, f{w,w) = 



fiw). From the definition of Radon transform, Eq. (19) 



f{w) = 27r / f{v)vdv. (21) 



We would have obtained the same result starting from Eq. (A12) 



B. Moving observer 

WIMP velocity distributions are often given in the galactic rest frame, while we are interested in the recoil momen- 
tum spectrum in the laboratory frame of the detector. The change of velocity frame can be performed either on the 
velocity distribution before computing the Radon transform or on the Radon transform computed in the galactic rest 
frame. The latter is often easier to compute, and the change of reference frame can be done simply as follows. 

The WIMP velocities viab and Vgai in the laboratory and galactic rest frames, respectively, are related by 

Vlab = Vgal - Viab, (22) 

where Viab is the velocity of the laboratory with respect to the galactic rest frame. This velocity transformation is a 



translation in velocity space, and we can use Eq. (A9) in Appendix A to relate the Radon transforms in the galactic 
and laboratory frames, 

/lab(w, w) = /gal (W + Vlab • W,w). (23) 

Thus the recoil momentum spectrum in the laboratory frame is given directly in terms of the Radon transform 
/gai(u', w) of the WIMP velocity distribution in the galactic rest frame by 

/gai(w9 + Viab-q, q), (24) 



dEdilq 47r^2 
with Vq — g/2/x as before. 

C. Rotated observer 

If we rotate the coordinate system, we see from Eq. ( [A8| ) in Appendix A that the recoil momentum spectrum is 
simply rotated, with the magnitude of the recoil momentum remaining the same, as expected. 
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V. EXAMPLES OF RECOIL MOMENTUM SPECTRA 

We give some examples of recoil momentum spectra corresponding to common velocity distributions. We obtain 
the recoil spectra for streams of particles and for isotropic and anisotropic Gaussian distributions with and without 
bulk velocities. 

A. A WIMP stream or flow 

The simplest case is that of a particle stream in which all WIMPs in the stream move with the same velocity V. 
In this case, 



and 



/stream (v) = (5(v - V), 



/stream (W,w) = (5( V • W - tu) . 



(25) 



(26) 



The recoil spectrum of a stream with velocity V is concentrated on a sphere of radius V jl, centered in V/2 and 
passing through the origin. The stream velocity V is a diameter of the sphere. 

Fig. ^ shows the {wx,Wy) section of the recoil momentum spectrum of a stream of WIMPs arriving from the left 
with velocity Vx = 400 kni/s. The full distribution is obtained through a rotation around the Wx axis. The pattern 
of recoil momenta forms a sphere. 



B. Maxwellian distribution 



A Maxwellian distribution with velocity dispersion , 



/m(w) = 



1 



■ exp 



V 

'2^ 



(27) 



(2^(72)3/2 

is a particular case of isotropic distribution, and we can use Eq. (pT|) above to compute its Radon transform. We find 



/m(w) 



1 



■ exp 



w 
'2^ 



(27ra2)i/2 

If the detector has velocity Viab, we can use Eq. (|23) to find the Radon transform in the laboratory frame. 



/M,lab(u', W) 



(2^2)1/2 



exp 



W ■ Vlab]' 

2a2 



(28) 



(29) 



Notice that w • Viab is the projection of the velocity of the observer in the direction of the nuclear recoil. This 
expression coincides with, but is simpler than, the analogous expression obtained by elementary methods in Ref. 

(cos 7 in Ref. ^ is cos 7 — — w • Vjab)- 

The recoil momentum distribution for a Maxwellian distribution is shown in Fig. |^, assuming a velocity dispersion 
of 300 km/s and an observer moving at 220 km/s in direction —x. The distribution is symmetric around the observer 
velocity. The figure shows the section in the {wx,Wy) plane only. The full distribution can be obtained by symmetry. 

Eq. ( p9| ) illustrates the reason for writing /(w, w) instead of /(w) (see Section A 1 for more details). The function 
/M.iab(0,w) assumes different values for different directions w; the function /(w) would therefore be multivalued at 
the origin. 



Truncated Maxwellian distribution 



We may truncate a Maxwellian distribution at the escape speed Vc 



1 



/tm(«) = { TVesc (271^2)3/2 

0, 



exp 



V 

2^. 



V < Wesc 

otherwise 



(30) 
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with 



erf 



/ Vesc \ 

V V2ay J 



2 



■ exp 



esc 



Then we have 



/tm(w) 



1 





r 2 1 




r 1 

CSC 




jexp 


.~2^. 


— exp 


."2a2_ 


} 



(31) 



(32) 



iVesc(27ra2)i/2 

D. Non-isotropic Gaussian distribution 

The recoil-momentum spectrum corresponding to an anisotropic Gaussian distribution can also be obtained ana- 
lytically. An anisotropic Gaussian distribution with variance matrix cr^ and mean velocity V is given by 



/gauss (v) 



1 



?7r3 detff^)!/ 



2u7?exp 



(v- V)^cr-2 (v- V) 



(33) 



We are using matrix notation, v"^ being the transpose of v, etc. Using the Fourier slice theorem, actually Eq. ( A14 ), 
we find the Radon transform of the anisotropic Gaussian to be 



/gauss(«^, W) 



1 



(27rw'rcr2 w)i/2 



exp 



2w'^(t2 w 



(34) 



This is another example of a function which assumes different values at w = according to the direction w. 



VI. RECONSTRUCTING THE VELOCITY DISTRIBUTION 



The recoil spectrum of a stream and a Maxwellian velocity distribution are very different: a sphere the first, a 
smooth distribution the second. This suggests that it may be possible to distinguish different kinds of WIMP velocity 
distributions just by examining the pattern of recoil momenta. Subtle differences among velocity distributions may 
be revealed by a maximum likelihood analysis of the corresponding recoil spectra |l^, |ll] . 

More ambitiously, we may think of recovering the WIMP velocity distribution by inverting the measured recoil 
momentum spectrum. Indeed, if we know the nuclear form factor of the detector nuclei, then for any fixed WIMP 
mass we can estimate the Radon transform of the WIMP velocity distribution from the measured recoil momentum 
spectrum, modulo a normalization constant K. Eq. ( |l7| ) can in fact be written as 

enabling us to obtain a measurement of the Radon transform f{vq,q) of the WIMP velocity distribution from the 
measured recoil spectrum dR/dEdflq. We may be able to invert this Radon transform and obtain the WIMP velocity 
distribution /(v), again modulo a normalization constant. Finally, we may be able to fix the normalization constant 
either by normalizing /(v) to unit integral or better by examining the detector efficiency as a function of WIMP 
velocity. 

There are several analytic formulas for the inversion of three-dimensional Radon transforms. Some of these formulas 
are collected in Appendix B for convenience. Most of the analytical inversion formulas can be converted into numerical 
algorithms |l9|] . However, any inversion algorithm we were able to find in the literature is suited only to a large amount 
of data in recoil momentum space, since they all assume that it would be possible to define a discretized version of 
/(t(;,w). This is not the case for directional dark matter searches, where the total number of events is not under the 
control of the experimentalist and is expected to be rather small. 

New inversion algorithms suited to small numbers of events are therefore needed if one wants to reconstruct the 
WIMP velocity distribution using data from directional detectors. As a first attempt in this direction, we have devised 
the following simple algorithm. Divide the WIMP velocity space into small cells Sm, to = 1, ... , M, and assume that 
the WIMP velocity distribution /(v) is constant over each of these small cells, with value in cell Sm- To each 
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recorded event j with measured recoil momentum q^, j = 1, . . . , iV, associate the plane Pj in WIMP velocity space 
defined by the equation 

Pj : 2fiv ■ clj = Qj. (36) 

Pj is the plane orthogonal to the recoil direction qj and at a distance Wj = qj/lji from the origin. Velocity vectors 
on this plane are all the WIMP velocities that can produce the observed nuclear recoil. Let 

= area(S'm n Pj), (37) 

the area of the intersection of the plane Pj with the cell Sm (see Appendix C for an explicit expression) . For each 

event j, assign weight ajm to the m-th cell. Sum the weights over the events, Am = X^jLi ^.j-m, essentially counting 
how many planes cross any given cell. Take the discrete Laplacian of the sum of the weights, and keep only those 
cells whose values exceed a predetermined threshold. The resulting distribution of cell values is our estimate of the 
WIMP velocity distribution. 

To test the capabilities of our algorithm, we simulated the recoil spectrum due to two streams of WIMPs arriving 
at the detector from opposite directions, with velocities {Vix,Viy,Viz) = (0,0,0.5) and {V2x,V2y,V2z) = (0,0,-0.2) 
(in arbitrary units). We generated 100 events, and applied the previous algorithm with 64^ cells in velocity space and 
a threshold of 0.1. We found that only two cells in velocity space are above threshold, and they correspond exactly 
to the location of the simulated streams. Fig. || plots the {vx,Vz) section of the reconstructed velocity distribution. 
It is impressive that we were able to recover this velocity distribution with only 100 events. 

We leave further studies of our simple algorithm, and the development of other algorithms, to future work. 



VII. CONCLUSIONS 



Directional detectors for WIMP dark matter searches will be able to measure not only the energy but also the 
direction of the nuclear recoils caused by the elastic scattering of galactic WIMPs with nuclei in the detector. This 
directional capability will help in separating a WIMP signal from background, and will also provide a measurement 
of the recoil momentum spectrum as compared to just the recoil energy spectrum. 

To gain insight into the properties of recoil momentum spectra, we have devised a simple formalism for the analytic 
computation of recoil momentum spectra from WIMP velocity distributions. Mathematically, the recoil momentum 
spectrum is the 3-dimensional Radon transform of the velocity distribution. 

Well-established mathematical properties of the Radon transform allow the computation of analytical expressions 
for recoil spectra associated to several common WIMP velocity distributions. As examples we presented recoil spectra 
for a WIMP stream, a Maxwellian, a truncated Maxwellian, and a non-isotropic Gaussian. We found in particular 
that a stream of WIMPs produces a characteristic spherical pattern of nuclear recoils. A Maxwellian distribution 
gives instead a smooth recoil pattern. Other velocity distributions lead to more complicated spectra. 

The analytic expressions we found for the nuclear recoil spectra will facilitate the discrimination of different velocity 
distributions through likelihood analysis. In addition, it may be possible to invert the measured momentum spectrum 
to reconstruct the local WIMP velocity distribution from data. For this purpose, we have presented an algorithm to 
recover the velocity distribution from a small number of recorded events. We have successfully recovered a simulated 
velocity distribution with just 100 generated events. 

We expect that the tools we have presented will be useful for the design and analysis of directional WIMP detectors. 
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APPENDIX A: SOME MATHEMATICS OF THE RADON TRANSFORM 



In this appendix we collect some useful mathematical properties of the 3-dimensional Radon transform. We denote 
the 3-dimensional Radon transform of a function /(v) by /(w,w). It is defined by 



f{w, w) = / (5(u; - w • w)f{v)d-'v. (Al) 
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It is easy to see that the Radon transform is Unear, 

f7+l2 = fl+f2- (A2) 
1. A remark on notation 

One may be tempted to write /(w) for /(w,w), after all w = ww. This notation may however be amb iguo us 



and should be used with care. Indeed, one must keep in mind that the Radon transform as defined in Eq. (Al) is 
a function of the magnitude w and the direction w separately. In other words, one may have /(0,w) 7^ /(0,w') for 
w 7^ w'. Namely, /(O, w) may assume different values for different directions. This will not be reflected in the notation 
/(w). The latter would read /(O) at the origin, independently of the direction w. In other words, /(w) would be a 
multiple-valued function at the origin. Mathematically, the distinction between /(w,w) and /(w) is important, and 
is expressed by saying that f{w,\v) is defined on M x S"^ while /(w) is defined on R'^. For our application, however, 
the distinction is of little concern, since the problematic origin w = corresponds to the region of vanishingly small 
recoil momenta, which is experimentally inaccessible. We have nevertheless used the mathematically correct notation 
throughout for clarity. 

2. Change of coordinates 

Under linear transformations of the coordinate axes, v transforms as 

V ^ v' = Av + b, (A3) 

where A is a 3 x 3 non-singular matrix and b is a constant vector. A velocity distribution function /(v) transforms 
so as to keep the number of particles in a volume (Pv invariant: 

f'iV)d'v' = fMd'v. (A4) 

Hence, 

/V) = ^/(A-(v'-b)), (A5) 
where det A is the determinant of A and A^^ is the inverse of A. To find the relation between the Radon transforms 



of /(v) and /'(v'), we change integration variable from v' to v in the definition, Eq. (Al), 

f'iw', w') = J 6{w' - w' • v') /'(v') d^v' = J 5{w' - w' • Av - w' • b) /(v) d^v 

= J S{w' - w • b - A^w • v) /(v) d^v = f{w' - w • b, A^w ), (A6) 
where A-^ is the transpose of A. Thus 

f'{w\ w') = f{w' - w' • b, A^w'). (A7) 

In particular, under a pure rotation R, 

f'iw',w')^f{w',R-'w'), (AS) 

and under a pure translation b, 

/'(u;',w') = .f(u''-w'-b, w'). (A9) 
3. Transformation of derivatives 



The following relations hold for derivatives of the Radon transform (here v ~ (vi, V2, W3) and w = (wi, W2, ^3)) 



4. Integration over angles 

We find the following expression for the integral of the Radon transform f{w, w) over the directions w: 



2tt j S{v cos^ — w) dcosj 



(5(v • w — w) dflu 



271 



/(v) d^v 



fiv)d-'v= I —eiv-\w\) f{v)d\ 



= 27r 



d^V. 



(A12) 



5. Fourier slice theorem 

There is a conne ctio n between the Radon transform and the Fourier transform. Taking the Fourier transform of 
the definition, Eq. (|Al[), with respect to w at fixed w gives 



+ 00 



rfAe*^7(A,w)= / f{v)e'''"-''d'v 



(A13) 



This equation goes imder the name of Fourier slice theorem. The right hand side is just the Fourier transform of /(v) 

evaluated at tw, while the left hand side is the Fourier transform of /(w, w) at fixed w. 
Inverting the Fourier transform in the left hand side of the Fourier slice theorem, we have 



f{w,w) = ^ / die-™* / /(v)e^*--d3^. 



(A14) 



This alternative expression of the Radon transform actually serves as its definition when functions are replaced by 
distributions (in the mathematical sense, see Ref. [EtI). 



6. Expansion into spherical harmonics 

Let us expand /(v) and its Radon transform /(w,w) into spherical harmonics Yijn{v) and y/rn(w), respectively. 
We have 



Im 



and 



The coefficients of the spherical harmonic expansions are related by 

fimiw) = 27r / fhn{v)vdv. 



(A15) 



(A16) 



(A17) 



where Pi{x) is a Legendre polynomial of order I. These expressions are useful when the velocity distributions are not 
isotrop ic. 

Eq. ( A17| ) can be proven using the decomposition of the (5-function in Legendre polynomials 



2/ + 1 

,5(v • w - <) = e{\ -t)Y, W Pii^ ■ w), 

the addition theorem for spherical harmonics 



(A18) 



(A19) 
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and the orthogonality of the spherical harmonics 

J Yir^,{^)Yimi^)da, = dl,lS.m'm, (A20) 

which lead to the relation 

' 5{v-w- t)Yira{^)d^v = 2tt9{1 - t)Pi{t)Yi,^{^). (A21) 



APPENDIX B: INVERSION FORMULAS FOR THE RADON TRANSFORM 
1. Inversion using the Laplacian 

An inversion formula for the Radon transform is 

/(^) = -i^ / /(v-w,w)dl]„, (Bl) 
where d^/dv'^ is the Laplacian in v. It can also be written in terms of /"(w,w) — f {w , w) / dw'^ as 



fi^) = -S^ J /"(v-w,w)dr!™, (B2) 

The inver sion formula ( |Bl[ ) can be proven by inverting the Fourier transform of /(v) in the Fourier slice theorem, 
Eq. ( [Al3| ), then integrating separately in t and w, and finally using the relation 



t^e''''dt^~TT6'^^\x), (B3) 
where (5^^^(x) is the second derivative of the Dirac J-function. 

2. Inversion through spherical harmonics 



Another inversion method is through an expansion in spherical harmonics. Referring to Eqs. (A_15) and ( [A16| ), one 
can prove the following inversion formula 

hn{v) = f Pi{-) fLHdw, (B4) 

where Pi{x) is a Legendre polynomial and f'il^{w) — d^fim/dw^, the second derivative of f{w) with respect to the 

modulus of w. 

Eq. (84) is proven along the same lines as Eq. ( |A17 ), starting from Eq. (B2) written as 



/(v) = / - V w)/"(u;, ^)dwdn^. (B5) 

3. Inversion through Fourier transforms 

The Fourier slice theorem, Eq. ( |A13| ), can be made into an algorithm for the numerical evaluation of the inverse 
Radon transform. Typically one would use fast Fourier transforms. 

4. Inversion through the adjoint operator 



Eq. (Bl) can also be made into an algorithm. For each given v, the integration in dil^ amounts to an integration over 
the sphere of diameter v and passing through the origin (a "stream sphere"), with integration measure dcos{9/2)d(p 
in spherical coordinates centered at the center of the sphere and north pole in v. The final Laplacian can be computed 
numerically as the difference between the central value and the average value of its six nearest neighbors. 
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5. Algebraic inversion via discretization 

An algebraic inversion method is the following JT?] , p. 96]. Suppose that the values fj,j — l,...,N, corresponding 
to the points Wj are known. In medical applications, the points Wj form a grid or other structure in space, and the 
/j's are the measured signal intensities. In our case, the number of detected events may be quite small, in which case 
we may let Wj be the actual measurement of a nuclear recoil momentum, with j varying over the number of events, 

and fj = l/sj, where Ej is the efficiency for detecting event j. 
By definition of Radon transform we have 

/ f{v)d'v = %, (B6) 

where the integral is taken over the plane in v-space defined by the equation 

Pj : v-w = w. (B7) 

Pj is the plane orthogonal to the recoil direction w and at a distance w from the origin. Now suppose that /(v) 
has compact support, meaning that it vanishes for |v| > something. This is a technical simplification that is valid in 
practice since real velocity distributions are always truncated at some large velocity (e.g. at the escape speed from 
the galaxy). Divide the v-space into small cells Sm, m = 1, . . . , M, and assume that /(v) is constant over each of 
these small cells, with value fm in cell Sm- This is the discretizing approximation. Let 

= area(S'™ n Pj), (B8) 

be the area of the intersection of the plane Pj with the cell Sm (see Appendix C for an explicit expression). A 



discretized version of Eq. ( B6 ) is then 



ajmfm. ^ fj- (B9) 



In matrix form 



Af = /, (BIO) 

where A = (ajm) is an iV x M matrix, / — (/i, . . . , JmY' / — (/ii ■ ■ • i In)'^- This is a system of linear equations 
for / that can be solved by inverting A. Since few ajm differ from zero, A is a sparse matrix, and it is convenient to 
solve this system iteratively. Fix uj, < uj < 2. Let the initial guess be /'^"^ and the fc-th update be /''^•'. From Z*^*^' 
compute the following vectors successively 

/(^) = f^'\ (Bll) 



/i'^ = /S + ^ [fi j «J ' J^h---,N. (B12) 

Here a| = J^m'^'jm ^^^^ '^j — ('^jii • ■ ■ i'^jm)- Finally let the next update be /C^+i) = Ref. |l^ attributes this 

method to Kaczmarz. 



UJ 



APPENDIX C: AREA OF THE INTERSECTION BETWEEN A PLANE AND A CELL 

For future reference, we give here the expression for the area of the intersection of a plane with a rectangular cell. 
Assume the {x,y,z) space is divided into rectangular cells of sides h^, hy, and along x, y, and z, respectively. 
Let the k)-th cell Sijk be centered in {xq + ih^^yo + jhy, zo + kh^), 

{i ~ \) hx + Xo < X < (i + 5) hx + Xq, 

S^Jk■. { {j -\)hy + yG<y <{i + \)hy + yQ, (CI) 
{k ~ \) + zq < z < (/c + i) ft,^ + zq. 



Let P be the plane defined by 



Lo^x + ujyy + uJzZ = p. (C2) 
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Then the area of the intersection of the plane P with the k)-th cell Sijk is 



area(5'yfe DP) — h^hyhz x < 



0, 
1 

1 P^+K^ 



J_ _ {K - P)^ 

[ kI ~ 2K1K2K3 ' 



if P,/K't + K'i + Kl > 73/2, 
if P < lifl and K>0, 

if P < l/Cl and K <0, 

ifP> 



(C3) 



where 



P = \p- uj^Xxq +ihx) - ujy{yo + jhy) - uj^{zo + khz)\, 
K = 1 (i^a - - i^i) , 



(C4) 
(C5) 



and -fCi, ^^2, and if 3 are t he q uantities l/ij^Wj^l, |/ij,ti;j,|, and |/izWz| sorted in order of increasing magnitude, Ki < K2 < 
K3. The last case in Eq. (|C^ ) becomes 



jf2 + iCs - 2P 

2if2if3 



(C6) 



in the limit of small Ki . 
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FIG. 1: Kinematics of elastic WIMP- nucleus scattering. 
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FIG. 2: Probability density distribution of the nuclear recoil momentum in the recoil plane (wxjWy), assuming a stream of 
WIMPs with velocity {vx,Vy,Vz) = (400km/s, 0, 0). The full {wx,Wy,Wz) distribution can be obtained by revolution around 
the Wx axis. The recoil momenta describe a sphere in recoil space. 
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FIG. 3: Probability density distribution of the nuclear recoil momentum in the recoil plane (wx,Wy), assuming a Maxwellian 
velocity distribution of WIMPs with velocity dispersion 300 km/s, and a detector moving with velocity {Vx,Vy,Vz) = 
(—220 km/s, 0, 0). Lighter areas have higher probability. The full (wxjWyjWz) distribution can be obtained by revolution 
around the Wx axis. 
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FIG. 4: Reconstructed velocity distribution of two WIMP streams with velocities (0, 0, 0.5) and (0, 0, —0.2) (in arbitrary units). 
Only the {vx,Vz) section is shown. 



